%% Q2: 把旧网格结点的数据插值到新网格

%% 读入数据
clear; clc;
timerVal = tic; % 对整个脚本计时
fileID = fopen('Grid_and_data.txt', 'r');
fmt = '%f %f %f %f';
A = textscan(fileID, fmt, 'HeaderLines', 1);
fclose(fileID);
S1 = A{4};

fileID = fopen('Result_grid_coordinates.txt', 'r');
fmt = '%d %f %f %f';
A = textscan(fileID, fmt, 'HeaderLines', 1);
fclose(fileID);
x = A{2}; y = A{3}; z = A{4};

% 读入数据计时
elapsed_time_readdata = toc(timerVal);
fprintf('读入数据用时：%f s\n', elapsed_time_readdata)

% 重新组织数据
Nx = 201; Ny = 101; Nz = 51;
% 生成网格向量 grid vectors
xg = linspace(0, 2000, Nx); yg = linspace(0, 1000, Ny); zg = linspace(0, 500, Nz);
% 将应力组织成三维网格的格式
S1_C=reshape(S1, [Ny, Nx, Nz]);

%% 结构网格插值
tic
F = griddedInterpolant({yg, xg, zg}, S1_C, 'linear');
S1_D = F(y, x, z);
elapsed_time = toc;
fprintf('插值用时：%f s\n', elapsed_time)

%% 输出到文件
tic
result = [double(A{1}) A{2} A{3} A{4} S1_D]';
fileID = fopen('Result_stress.txt', 'w');
fprintf(fileID, 'No.\tx\ty\tz\tS1\r\n');
fprintf(fileID, '%8d %6.2f %6.2f %6.2f %12.8f\r\n', result);
fclose(fileID);
elapsed_time = toc;
fprintf('输出用时：%f s\n', elapsed_time)

%% 总用时
elapsed_time2 = toc(timerVal);
fprintf('总用时：%f s\n', elapsed_time2)